//
// Created by Lenovo on 2024/11/13.
//



#include <complex>
#include <iostream>
#include <cmath>
#include <vector>

constexpr  int MAX_N = 1 << 20;
constexpr double PI = 3.141259;

using Complex = std::complex<double>;
std::complex<double> tmp[MAX_N];

void DFT(std::complex<double> *f, int n, int rev)
{
    if (1 == n)
        return ;
    double PI = acos(-1);

    using Comp = std::complex<double>;
    for (int i = 0;i < n;++i)
        tmp[i] = f[i];
    for (int i = 0;i < n;++i) {
        if ( i & 1)
            f[n/2 + i / 2] = tmp[i];
        else
            f[i/2] = tmp[i];
    }
    Comp *e = f;
    Comp *o = f + n / 2;

    DFT(e, n / 2, rev);
    DFT(o, n / 2, rev);

    Comp cur(1, 0);
    Comp step(std::cos(2 * PI/n), std::sin(2 * PI *rev/n));
    for (int k = 0;k < n / 2; ++k) {
        tmp[k] = e[k] + cur * o[k];
        tmp[k + n / 2] = e[k] - cur * o[k];
        cur *= step;
    }
    for (int i = 0;i < n; ++i)
            f[i] = tmp[i];
}

// 同样需要保证 len 是 2 的幂
// 记 rev[i] 为 i 翻转后的值
void change(Complex y[], int len) {

    std::vector<int> rev(MAX_N, 0);
    for (int i = 0; i < len; ++i) {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1) {  // 如果最后一位是 1，则翻转成 len/2
            rev[i] |= len >> 1;
        }
    }
    for (int i = 0; i < len; ++i) {
        if (i < rev[i]) {  // 保证每对数只翻转一次
            swap(y[i], y[rev[i]]);
        }
    }
    return;
}

// Ref: OI-wiki-fft!
void fft(Complex y[], int len, int on) {
    // 位逆序置换
    change(y, len);
    // 模拟合并过程，一开始，从长度为一合并到长度为二，一直合并到长度为 len。
    for (int h = 2; h <= len; h <<= 1) {
        // wn：当前单位复根的间隔：w^1_h
        Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
        // 合并，共 len / h 次。
        for (int j = 0; j < len; j += h) {
            // 计算当前单位复根，一开始是 1 = w^0_n，之后是以 wn 为间隔递增： w^1_n
            // ...
            Complex w(1, 0);
            for (int k = j; k < j + h / 2; k++) {
                // 左侧部分和右侧是子问题的解
                Complex u = y[k];
                Complex t = w * y[k + h / 2];
                // 这就是把两部分分治的结果加起来
                y[k] = u + t;
                y[k + h / 2] = u - t;
                // 后半个 「step」 中的ω一定和 「前半个」 中的成相反数
                // 「红圈」上的点转一整圈「转回来」，转半圈正好转成相反数
                // 一个数相反数的平方与这个数自身的平方相等
                w = w * wn;
            }
        }
    }
    // 如果是 IDFT，它的逆矩阵的每一个元素不只是原元素取倒数，还要除以长度 len。
    if (on == -1) {
        for (int i = 0; i < len; i++) {
            y[i] /= len;
        }
    }
}



